!-------------------------------------------------------------LICENSE--------------------------------------------------------------!
!                                                                                                                                  !
!The MAP code is written in Fortran language for magnetohydrodynamics (MHD) calculation with the adaptive mesh refinement (AMR)    !
!and Message Passing Interface (MPI) parallelization.                                                                              !
!                                                                                                                                  !
!Copyright (C) 2012                                                                                                                !
!Ronglin Jiang                                                                                                                     !
!rljiang@ssc.net.cn                                                                                                                !
!585 Guoshoujing Road. Pudong, Shanghai, P.R.C. 201203                                                                             !
!                                                                                                                                  !
!This program is free software; you can redistribute it and/or modify it under the terms of the GNU General Public License         !
!as published by the Free Software Foundation; either version 2 of the License, or (at your option) any later version.             !
!                                                                                                                                  !
!This program is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the implied warranty of    !
!MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License for more details.                        !
!                                                                                                                                  !
!You should have received a copy of the GNU General Public License along with this program; if not, write to the Free Software     !
!Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301, USA.                                                   !
!                                                                                                                                  !
!-------------------------------------------------------------LICENSE--------------------------------------------------------------!

!==================================================================================================================================|
subroutine block_reflux (current)
!==================================================================================================================================|

   use parameters
   use block_define
   use interface_block_boundary_physi
   implicit none

   type(block), pointer, intent(in) :: current

   integer(4) :: i_n, j_n, k_n, i_block, j_block, k_block, i_neigh, j_neigh, k_neigh, i_neigh_n, j_neigh_n, k_neigh_n
   integer(4) :: i2, j2, k2, i2m1, j2m1, k2m1, n_reflux
   integer(4) :: nx_tmp, ny_tmp, nz_tmp, i_end, j_end, k_end
   real(8) :: dt_dx, dt_dy, dt_dz

   type(block), pointer :: neigh, child

!----------------------------------------------------------------------------------------------------------------------------------|
   dt_dx = min_dt / current%dx
   dt_dy = min_dt / current%dy
   dt_dz = min_dt / current%dz

!----------------------------------------------------------------------------------------------------------------------------------|
!  left boundary
!----------------------------------------------------------------------------------------------------------------------------------|
   i_neigh = -1
   j_neigh = 0
   k_neigh = 0
   neigh => current%neigh(i_neigh, j_neigh, k_neigh)%pblock
   if (associated(neigh)) then
      if (current%nc .eq. 0 .and. neigh%nc .ne. 0 .and. reflux_flag .eq. 1) then
         i_block = neigh%nc_x
         i_neigh = 1
         i_neigh_n = 2
         do k_block = 1, neigh%nc_z
            k_end = nz_no_bnd / ratio
            nz_tmp = (k_block - 1) * nz_no_bnd / ratio + ng
            do j_block = 1, neigh%nc_y
               j_end = ny_no_bnd / ratio
               ny_tmp = (j_block - 1) * ny_no_bnd / ratio + ng
               child => neigh%child(i_block, j_block, k_block)%pblock
               if (child%cpu_id .ne. cpu_id) cycle
               n_reflux = ng + 1
               do k = 1, k_end
                  k_n = k + nz_tmp
                  k2 = k * ratio + ng
                  k2m1 = (k - 1) * ratio + 1 + ng
                  do j = 1, j_end
                     j_n = j + ny_tmp
                     j2 = j * ratio + ng
                     j2m1 = (j - 1) * ratio + 1 + ng
                     current%ro(n_reflux, j_n, k_n) = current%ro(n_reflux, j_n, k_n) + dt_dx *                                     &
                        (sum(child%fxro_bnd(j2m1 : j2, k2m1 : k2, i_neigh_n)) / ratio2 - current%fxro_bnd(j_n, k_n, i_neigh))
                     current%mx(n_reflux, j_n, k_n) = current%mx(n_reflux, j_n, k_n) + dt_dx *                                     &
                        (sum(child%fxmx_bnd(j2m1 : j2, k2m1 : k2, i_neigh_n)) / ratio2 - current%fxmx_bnd(j_n, k_n, i_neigh))
                     current%my(n_reflux, j_n, k_n) = current%my(n_reflux, j_n, k_n) + dt_dx *                                     &
                        (sum(child%fxmy_bnd(j2m1 : j2, k2m1 : k2, i_neigh_n)) / ratio2 - current%fxmy_bnd(j_n, k_n, i_neigh))
                     current%mz(n_reflux, j_n, k_n) = current%mz(n_reflux, j_n, k_n) + dt_dx *                                     &
                        (sum(child%fxmz_bnd(j2m1 : j2, k2m1 : k2, i_neigh_n)) / ratio2 - current%fxmz_bnd(j_n, k_n, i_neigh))
                     current%bx(n_reflux, j_n, k_n) = current%bx(n_reflux, j_n, k_n) + dt_dx *                                     &
                        (sum(child%fxbx_bnd(j2m1 : j2, k2m1 : k2, i_neigh_n)) / ratio2 - current%fxbx_bnd(j_n, k_n, i_neigh))
                     current%by(n_reflux, j_n, k_n) = current%by(n_reflux, j_n, k_n) + dt_dx *                                     &
                        (sum(child%fxby_bnd(j2m1 : j2, k2m1 : k2, i_neigh_n)) / ratio2 - current%fxby_bnd(j_n, k_n, i_neigh))
                     current%bz(n_reflux, j_n, k_n) = current%bz(n_reflux, j_n, k_n) + dt_dx *                                     &
                        (sum(child%fxbz_bnd(j2m1 : j2, k2m1 : k2, i_neigh_n)) / ratio2 - current%fxbz_bnd(j_n, k_n, i_neigh))
                     current%en(n_reflux, j_n, k_n) = current%en(n_reflux, j_n, k_n) + dt_dx *                                     &
                        (sum(child%fxen_bnd(j2m1 : j2, k2m1 : k2, i_neigh_n)) / ratio2 - current%fxen_bnd(j_n, k_n, i_neigh))
      !               current%po(n_reflux, j_n) = current%po(n_reflux, j_n) + dt_dx *                                               &
      !                  (sum(child%fxpo_bnd(j2m1 : j2, i_neigh_n)) / ratio - current%fxpo_bnd(j_n, i_neigh))
                  enddo
               enddo
            enddo
         enddo
      endif
   endif

!----------------------------------------------------------------------------------------------------------------------------------|
!  right boundary
!----------------------------------------------------------------------------------------------------------------------------------|
   i_neigh = 1
   j_neigh = 0
   k_neigh = 0
   neigh => current%neigh(i_neigh, j_neigh, k_neigh)%pblock
   if (associated(neigh)) then
      if (current%nc .eq. 0 .and. neigh%nc .ne. 0 .and. reflux_flag .eq. 1) then
         i_block = 1
         i_neigh = 2
         i_neigh_n = 1
         do k_block = 1, neigh%nc_z
            k_end = nz_no_bnd / ratio
            nz_tmp = (k_block - 1) * nz_no_bnd / ratio + ng
            do j_block = 1, neigh%nc_y
               j_end = ny_no_bnd / ratio
               ny_tmp = (j_block - 1) * ny_no_bnd / ratio + ng
               child => neigh%child(i_block, j_block, k_block)%pblock
               if (child%cpu_id .ne. cpu_id) cycle
               n_reflux = nx_block - ng
               do k = 1, k_end
                  k_n = k + nz_tmp
                  k2 = k * ratio + ng
                  k2m1 = (k - 1) * ratio + 1 + ng
                  do j = 1, j_end
                     j_n = j + ny_tmp
                     j2 = j * ratio + ng
                     j2m1 = (j - 1) * ratio + 1 + ng
                     current%ro(n_reflux, j_n, k_n) = current%ro(n_reflux, j_n, k_n) - dt_dx *                                     &
                        (sum(child%fxro_bnd(j2m1 : j2, k2m1 : k2, i_neigh_n)) / ratio2 - current%fxro_bnd(j_n, k_n, i_neigh))
                     current%mx(n_reflux, j_n, k_n) = current%mx(n_reflux, j_n, k_n) - dt_dx *                                     &
                        (sum(child%fxmx_bnd(j2m1 : j2, k2m1 : k2, i_neigh_n)) / ratio2 - current%fxmx_bnd(j_n, k_n, i_neigh))
                     current%my(n_reflux, j_n, k_n) = current%my(n_reflux, j_n, k_n) - dt_dx *                                     &
                        (sum(child%fxmy_bnd(j2m1 : j2, k2m1 : k2, i_neigh_n)) / ratio2 - current%fxmy_bnd(j_n, k_n, i_neigh))
                     current%mz(n_reflux, j_n, k_n) = current%mz(n_reflux, j_n, k_n) - dt_dx *                                     &
                        (sum(child%fxmz_bnd(j2m1 : j2, k2m1 : k2, i_neigh_n)) / ratio2 - current%fxmz_bnd(j_n, k_n, i_neigh))
                     current%bx(n_reflux, j_n, k_n) = current%bx(n_reflux, j_n, k_n) - dt_dx *                                     &
                        (sum(child%fxbx_bnd(j2m1 : j2, k2m1 : k2, i_neigh_n)) / ratio2 - current%fxbx_bnd(j_n, k_n, i_neigh))
                     current%by(n_reflux, j_n, k_n) = current%by(n_reflux, j_n, k_n) - dt_dx *                                     &
                        (sum(child%fxby_bnd(j2m1 : j2, k2m1 : k2, i_neigh_n)) / ratio2 - current%fxby_bnd(j_n, k_n, i_neigh))
                     current%bz(n_reflux, j_n, k_n) = current%bz(n_reflux, j_n, k_n) - dt_dx *                                     &
                        (sum(child%fxbz_bnd(j2m1 : j2, k2m1 : k2, i_neigh_n)) / ratio2 - current%fxbz_bnd(j_n, k_n, i_neigh))
                     current%en(n_reflux, j_n, k_n) = current%en(n_reflux, j_n, k_n) - dt_dx *                                     &
                        (sum(child%fxen_bnd(j2m1 : j2, k2m1 : k2, i_neigh_n)) / ratio2 - current%fxen_bnd(j_n, k_n, i_neigh))
      !               current%po(n_reflux, j_n) = current%po(n_reflux, j_n) - dt_dx *                                               &
      !                  (sum(child%fxpo_bnd(j2m1 : j2, i_neigh_n)) / ratio - current%fxpo_bnd(j_n, i_neigh))
                  enddo
               enddo
            enddo
         enddo
      endif
   endif

!----------------------------------------------------------------------------------------------------------------------------------|
!  behind boundary
!----------------------------------------------------------------------------------------------------------------------------------|
   i_neigh = 0
   j_neigh = -1
   k_neigh = 0
   neigh => current%neigh(i_neigh, j_neigh, k_neigh)%pblock
   if (associated(neigh)) then
      if (current%nc .eq. 0 .and. neigh%nc .ne. 0 .and. reflux_flag .eq. 1) then
         j_block = neigh%nc_y
         j_neigh = 1
         j_neigh_n = 2
         do k_block = 1, neigh%nc_z
            k_end = nz_no_bnd / ratio
            nz_tmp = (k_block - 1) * nz_no_bnd / ratio + ng
            do i_block = 1, neigh%nc_x
               i_end = nx_no_bnd / ratio
               nx_tmp = (i_block - 1) * nx_no_bnd / ratio + ng
               child => neigh%child(i_block, j_block, k_block)%pblock
               if (child%cpu_id .ne. cpu_id) cycle
               n_reflux = ng + 1
               do k = 1, k_end
                  k_n = k + nz_tmp
                  k2 = k * ratio + ng
                  k2m1 = (k - 1) * ratio + 1 + ng
                  do i = 1, i_end
                     i_n = i + nx_tmp
                     i2 = i * ratio + ng
                     i2m1 = (i - 1) * ratio + 1 + ng
                     current%ro(i_n, n_reflux, k_n) = current%ro(i_n, n_reflux, k_n) + dt_dy *                                     &
                        (sum(child%fyro_bnd(i2m1 : i2, k2m1 : k2, j_neigh_n)) / ratio2 - current%fyro_bnd(i_n, k_n, j_neigh))
                     current%mx(i_n, n_reflux, k_n) = current%mx(i_n, n_reflux, k_n) + dt_dy *                                     &
                        (sum(child%fymx_bnd(i2m1 : i2, k2m1 : k2, j_neigh_n)) / ratio2 - current%fymx_bnd(i_n, k_n, j_neigh))
                     current%my(i_n, n_reflux, k_n) = current%my(i_n, n_reflux, k_n) + dt_dy *                                     &
                        (sum(child%fymy_bnd(i2m1 : i2, k2m1 : k2, j_neigh_n)) / ratio2 - current%fymy_bnd(i_n, k_n, j_neigh))
                     current%mz(i_n, n_reflux, k_n) = current%mz(i_n, n_reflux, k_n) + dt_dy *                                     &
                        (sum(child%fymz_bnd(i2m1 : i2, k2m1 : k2, j_neigh_n)) / ratio2 - current%fymz_bnd(i_n, k_n, j_neigh))
                     current%bx(i_n, n_reflux, k_n) = current%bx(i_n, n_reflux, k_n) + dt_dy *                                     &
                        (sum(child%fybx_bnd(i2m1 : i2, k2m1 : k2, j_neigh_n)) / ratio2 - current%fybx_bnd(i_n, k_n, j_neigh))
                     current%by(i_n, n_reflux, k_n) = current%by(i_n, n_reflux, k_n) + dt_dy *                                     &
                        (sum(child%fyby_bnd(i2m1 : i2, k2m1 : k2, j_neigh_n)) / ratio2 - current%fyby_bnd(i_n, k_n, j_neigh))
                     current%bz(i_n, n_reflux, k_n) = current%bz(i_n, n_reflux, k_n) + dt_dy *                                     &
                        (sum(child%fybz_bnd(i2m1 : i2, k2m1 : k2, j_neigh_n)) / ratio2 - current%fybz_bnd(i_n, k_n, j_neigh))
                     current%en(i_n, n_reflux, k_n) = current%en(i_n, n_reflux, k_n) + dt_dy *                                     &
                        (sum(child%fyen_bnd(i2m1 : i2, k2m1 : k2, j_neigh_n)) / ratio2 - current%fyen_bnd(i_n, k_n, j_neigh))
      !               current%po(i_n, n_reflux) = current%po(i_n, n_reflux) + dt_dy *                                               &
      !                  (sum(child%fypo_bnd(i2m1 : i2, j_neigh_n)) / ratio - current%fypo_bnd(i_n, j_neigh))
                  enddo
               enddo
            enddo
         enddo
      endif
   endif

!----------------------------------------------------------------------------------------------------------------------------------|
!  front boundary
!----------------------------------------------------------------------------------------------------------------------------------|
   i_neigh = 0
   j_neigh = 1
   k_neigh = 0
   neigh => current%neigh(i_neigh, j_neigh, k_neigh)%pblock
   if (associated(neigh)) then
      if (current%nc .eq. 0 .and. neigh%nc .ne. 0 .and. reflux_flag .eq. 1) then
         j_block = 1
         j_neigh = 2
         j_neigh_n = 1
         do k_block = 1, neigh%nc_z
            k_end = nz_no_bnd / ratio
            nz_tmp = (k_block - 1) * nz_no_bnd / ratio + ng
            do i_block = 1, neigh%nc_x
               i_end = nx_no_bnd / ratio
               nx_tmp = (i_block - 1) * nx_no_bnd / ratio + ng
               child => neigh%child(i_block, j_block, k_block)%pblock
               if (child%cpu_id .ne. cpu_id) cycle
               n_reflux = ny_block - ng
               do k = 1, k_end
                  k_n = k + nz_tmp
                  k2 = k * ratio + ng
                  k2m1 = (k - 1) * ratio + 1 + ng
                  do i = 1, i_end
                     i_n = i + nx_tmp
                     i2 = i * ratio + ng
                     i2m1 = (i - 1) * ratio + 1 + ng
                     current%ro(i_n, n_reflux, k_n) = current%ro(i_n, n_reflux, k_n) - dt_dy *                                     &
                        (sum(child%fyro_bnd(i2m1 : i2, k2m1 : k2, j_neigh_n)) / ratio2 - current%fyro_bnd(i_n, k_n, j_neigh))
                     current%mx(i_n, n_reflux, k_n) = current%mx(i_n, n_reflux, k_n) - dt_dy *                                     &
                        (sum(child%fymx_bnd(i2m1 : i2, k2m1 : k2, j_neigh_n)) / ratio2 - current%fymx_bnd(i_n, k_n, j_neigh))
                     current%my(i_n, n_reflux, k_n) = current%my(i_n, n_reflux, k_n) - dt_dy *                                     &
                        (sum(child%fymy_bnd(i2m1 : i2, k2m1 : k2, j_neigh_n)) / ratio2 - current%fymy_bnd(i_n, k_n, j_neigh))
                     current%mz(i_n, n_reflux, k_n) = current%mz(i_n, n_reflux, k_n) - dt_dy *                                     &
                        (sum(child%fymz_bnd(i2m1 : i2, k2m1 : k2, j_neigh_n)) / ratio2 - current%fymz_bnd(i_n, k_n, j_neigh))
                     current%bx(i_n, n_reflux, k_n) = current%bx(i_n, n_reflux, k_n) - dt_dy *                                     &
                        (sum(child%fybx_bnd(i2m1 : i2, k2m1 : k2, j_neigh_n)) / ratio2 - current%fybx_bnd(i_n, k_n, j_neigh))
                     current%by(i_n, n_reflux, k_n) = current%by(i_n, n_reflux, k_n) - dt_dy *                                     &
                        (sum(child%fyby_bnd(i2m1 : i2, k2m1 : k2, j_neigh_n)) / ratio2 - current%fyby_bnd(i_n, k_n, j_neigh))
                     current%bz(i_n, n_reflux, k_n) = current%bz(i_n, n_reflux, k_n) - dt_dy *                                     &
                        (sum(child%fybz_bnd(i2m1 : i2, k2m1 : k2, j_neigh_n)) / ratio2 - current%fybz_bnd(i_n, k_n, j_neigh))
                     current%en(i_n, n_reflux, k_n) = current%en(i_n, n_reflux, k_n) - dt_dy *                                     &
                        (sum(child%fyen_bnd(i2m1 : i2, k2m1 : k2, j_neigh_n)) / ratio2 - current%fyen_bnd(i_n, k_n, j_neigh))
      !               current%po(i_n, n_reflux) = current%po(i_n, n_reflux) - dt_dy *                                               &
      !                  (sum(child%fypo_bnd(i2m1 : i2, j_neigh_n)) / ratio - current%fypo_bnd(i_n, j_neigh))
                  enddo
               enddo
            enddo
         enddo
      endif
   endif

!----------------------------------------------------------------------------------------------------------------------------------|
!  bottom boundary
!----------------------------------------------------------------------------------------------------------------------------------|
   i_neigh = 0
   j_neigh = 0
   k_neigh = -1
   neigh => current%neigh(i_neigh, j_neigh, k_neigh)%pblock
   if (associated(neigh)) then
      if (current%nc .eq. 0 .and. neigh%nc .ne. 0 .and. reflux_flag .eq. 1) then
         k_block = neigh%nc_z
         k_neigh = 1
         k_neigh_n = 2
         do j_block = 1, neigh%nc_y
            j_end = ny_no_bnd / ratio
            ny_tmp = (j_block - 1) * ny_no_bnd / ratio + ng
            do i_block = 1, neigh%nc_x
               i_end = nx_no_bnd / ratio
               nx_tmp = (i_block - 1) * nx_no_bnd / ratio + ng
               child => neigh%child(i_block, j_block, k_block)%pblock
               if (child%cpu_id .ne. cpu_id) cycle
               n_reflux = ng + 1
               do j = 1, j_end
                  j_n = j + ny_tmp
                  j2 = j * ratio + ng
                  j2m1 = (j - 1) * ratio + 1 + ng
                  do i = 1, i_end
                     i_n = i + nx_tmp
                     i2 = i * ratio + ng
                     i2m1 = (i - 1) * ratio + 1 + ng
                     current%ro(i_n, j_n, n_reflux) = current%ro(i_n, j_n, n_reflux) + dt_dz *                                     &
                        (sum(child%fzro_bnd(i2m1 : i2, j2m1 : j2, k_neigh_n)) / ratio2 - current%fzro_bnd(i_n, j_n, k_neigh))
                     current%mx(i_n, j_n, n_reflux) = current%mx(i_n, j_n, n_reflux) + dt_dz *                                     &
                        (sum(child%fzmx_bnd(i2m1 : i2, j2m1 : j2, k_neigh_n)) / ratio2 - current%fzmx_bnd(i_n, j_n, k_neigh))
                     current%my(i_n, j_n, n_reflux) = current%my(i_n, j_n, n_reflux) + dt_dz *                                     &
                        (sum(child%fzmy_bnd(i2m1 : i2, j2m1 : j2, k_neigh_n)) / ratio2 - current%fzmy_bnd(i_n, j_n, k_neigh))
                     current%mz(i_n, j_n, n_reflux) = current%mz(i_n, j_n, n_reflux) + dt_dz *                                     &
                        (sum(child%fzmz_bnd(i2m1 : i2, j2m1 : j2, k_neigh_n)) / ratio2 - current%fzmz_bnd(i_n, j_n, k_neigh))
                     current%bx(i_n, j_n, n_reflux) = current%bx(i_n, j_n, n_reflux) + dt_dz *                                     &
                        (sum(child%fzbx_bnd(i2m1 : i2, j2m1 : j2, k_neigh_n)) / ratio2 - current%fzbx_bnd(i_n, j_n, k_neigh))
                     current%by(i_n, j_n, n_reflux) = current%by(i_n, j_n, n_reflux) + dt_dz *                                     &
                        (sum(child%fzby_bnd(i2m1 : i2, j2m1 : j2, k_neigh_n)) / ratio2 - current%fzby_bnd(i_n, j_n, k_neigh))
                     current%bz(i_n, j_n, n_reflux) = current%bz(i_n, j_n, n_reflux) + dt_dz *                                     &
                        (sum(child%fzbz_bnd(i2m1 : i2, j2m1 : j2, k_neigh_n)) / ratio2 - current%fzbz_bnd(i_n, j_n, k_neigh))
                     current%en(i_n, j_n, n_reflux) = current%en(i_n, j_n, n_reflux) + dt_dz *                                     &
                        (sum(child%fzen_bnd(i2m1 : i2, j2m1 : j2, k_neigh_n)) / ratio2 - current%fzen_bnd(i_n, j_n, k_neigh))
      !               current%po(i_n, n_reflux) = current%po(i_n, n_reflux) + dt_dz *                                               &
      !                  (sum(child%fzpo_bnd(i2m1 : i2, j_neigh_n)) / ratio - current%fzpo_bnd(i_n, j_neigh))
                  enddo
               enddo
            enddo
         enddo
      endif
   endif

!----------------------------------------------------------------------------------------------------------------------------------|
!  top boundary
!----------------------------------------------------------------------------------------------------------------------------------|
   i_neigh = 0
   j_neigh = 0
   k_neigh = 1
   neigh => current%neigh(i_neigh, j_neigh, k_neigh)%pblock
   if (associated(neigh)) then
      if (current%nc .eq. 0 .and. neigh%nc .ne. 0 .and. reflux_flag .eq. 1) then
         k_block = 1
         k_neigh = 2
         k_neigh_n = 1
         do j_block = 1, neigh%nc_y
            j_end = ny_no_bnd / ratio
            ny_tmp = (j_block - 1) * ny_no_bnd / ratio + ng
            do i_block = 1, neigh%nc_x
               i_end = nx_no_bnd / ratio
               nx_tmp = (i_block - 1) * nx_no_bnd / ratio + ng
               child => neigh%child(i_block, j_block, k_block)%pblock
               if (child%cpu_id .ne. cpu_id) cycle
               n_reflux = nz_block - ng
               do j = 1, j_end
                  j_n = j + ny_tmp
                  j2 = j * ratio + ng
                  j2m1 = (j - 1) * ratio + 1 + ng
                  do i = 1, i_end
                     i_n = i + nx_tmp
                     i2 = i * ratio + ng
                     i2m1 = (i - 1) * ratio + 1 + ng
                     current%ro(i_n, j_n, n_reflux) = current%ro(i_n, j_n, n_reflux) - dt_dz *                                     &
                        (sum(child%fzro_bnd(i2m1 : i2, j2m1 : j2, k_neigh_n)) / ratio2 - current%fzro_bnd(i_n, j_n, k_neigh))
                     current%mx(i_n, j_n, n_reflux) = current%mx(i_n, j_n, n_reflux) - dt_dz *                                     &
                        (sum(child%fzmx_bnd(i2m1 : i2, j2m1 : j2, k_neigh_n)) / ratio2 - current%fzmx_bnd(i_n, j_n, k_neigh))
                     current%my(i_n, j_n, n_reflux) = current%my(i_n, j_n, n_reflux) - dt_dz *                                     &
                        (sum(child%fzmy_bnd(i2m1 : i2, j2m1 : j2, k_neigh_n)) / ratio2 - current%fzmy_bnd(i_n, j_n, k_neigh))
                     current%mz(i_n, j_n, n_reflux) = current%mz(i_n, j_n, n_reflux) - dt_dz *                                     &
                        (sum(child%fzmz_bnd(i2m1 : i2, j2m1 : j2, k_neigh_n)) / ratio2 - current%fzmz_bnd(i_n, j_n, k_neigh))
                     current%bx(i_n, j_n, n_reflux) = current%bx(i_n, j_n, n_reflux) - dt_dz *                                     &
                        (sum(child%fzbx_bnd(i2m1 : i2, j2m1 : j2, k_neigh_n)) / ratio2 - current%fzbx_bnd(i_n, j_n, k_neigh))
                     current%by(i_n, j_n, n_reflux) = current%by(i_n, j_n, n_reflux) - dt_dz *                                     &
                        (sum(child%fzby_bnd(i2m1 : i2, j2m1 : j2, k_neigh_n)) / ratio2 - current%fzby_bnd(i_n, j_n, k_neigh))
                     current%bz(i_n, j_n, n_reflux) = current%bz(i_n, j_n, n_reflux) - dt_dz *                                     &
                        (sum(child%fzbz_bnd(i2m1 : i2, j2m1 : j2, k_neigh_n)) / ratio2 - current%fzbz_bnd(i_n, j_n, k_neigh))
                     current%en(i_n, j_n, n_reflux) = current%en(i_n, j_n, n_reflux) - dt_dz *                                     &
                        (sum(child%fzen_bnd(i2m1 : i2, j2m1 : j2, k_neigh_n)) / ratio2 - current%fzen_bnd(i_n, j_n, k_neigh))
      !               current%po(i_n, n_reflux) = current%po(i_n, n_reflux) - dt_dz *                                               &
      !                  (sum(child%fzpo_bnd(i2m1 : i2, k_neigh_n)) / ratio - current%fzpo_bnd(i_n, k_neigh))
                  enddo
               enddo
            enddo
         enddo
      endif
   endif

!----------------------------------------------------------------------------------------------------------------------------------|
   return
end subroutine block_reflux
